#####################################################################
# Figure 2
#####################################################################

#####################################################################
# Preliminaries
#####################################################################

rm(list = ls())

library(sandwich)

set.seed(12345)

#####################################################################
# Functions
#####################################################################

source("code/auxiliary scripts/common_functions.R")

sim_nw <- function(T0,T1,rho,K,sig_level){

  # Setup
  
  T01   <- T0 + T1
  u     <- generate.AR.series(T01,rho)
  
  r       <- min(floor(T0/K),T1)
  u.pre   <- u[1:T0]
  u.post  <- u[(T0+1):T01]

  # t-Test
  
  tau.mat <- matrix(NA,K,1)
  for (k in 1:K){
    Hk            <- (T0-(r*K))+seq((k-1)*r+1,k*r,1) 
    tau.mat[k,1]  <- mean(u.post) - mean(u.pre[Hk])
  }
  
  tau.hat <- mean(tau.mat)
  se.hat  <- sqrt(1+((K*r)/T1))*sd(tau.mat)/sqrt(K)
  t.hat   <- tau.hat/se.hat # this has a t_{K-1} distribution
  
  cov.t   <- (abs(t.hat)<= qt(1-sig_level/2,df=K-1))
  
  # Newey-West
  
  m       <- lm(u.pre~ 1)
  lrv     <- NeweyWest(m, prewhite = TRUE)
  c0      <- T0/T1
  g.c0.K  <- K*(c0<1)+K/c0*(1<=c0&c0<=K)+1*(c0>K)
  std.lrv <- sqrt(lrv*T0*(min(c0,1)+g.c0.K/K))/sqrt(min(T0,T1))
  z.hat   <- tau.hat/std.lrv
  cov.z   <- (abs(z.hat)<= qnorm(1-sig_level/2))
  results <- c(cov.t,cov.z)
  
  return(results)
  
}

#####################################################################
# Simulations
#####################################################################

# Setup

nreps <- 50000
rho.vec <- c(seq(0,0.9,0.1))

# Simulation loop 
# Comment: The NeweyWest() package will generate the following warning for a few draws 
# of the data due to the very small sample size,
# "In meatHAC(x, order.by = order.by, prewhite = prewhite, weights = weights,  :
# more weights than observations, only first n used", 
# indicating automatic truncation of vector of weights to match the number of observations.
# The code before and after the simulation loop suppresses these warnings.

old_warning_option <- getOption("warn") 
options(warn = -1)

results_overall_t_15  <- results_overall_nw_15 <- results_overall_t_30  <- results_overall_nw_30 <- matrix(NA,length(rho.vec),1)
for (r in 1:length(rho.vec)){
  temp_mat_15 <- temp_mat_30 <- matrix(NA,nreps,2)
  for (rep in 1:nreps){
    temp_mat_15[rep,] <- sim_nw(T0=15,T1=16,rho=rho.vec[r],K=3,sig_level=0.1)
    temp_mat_30[rep,] <- sim_nw(T0=30,T1=16,rho=rho.vec[r],K=3,sig_level=0.1)
  }
  results_overall_t_15[r,]    <- mean(temp_mat_15[,1])
  results_overall_nw_15[r,]   <- mean(temp_mat_15[,2])
  results_overall_t_30[r,]    <- mean(temp_mat_30[,1])
  results_overall_nw_30[r,]   <- mean(temp_mat_30[,2])
}

options(warn=old_warning_option) 

#####################################################################
# Graphics
#####################################################################

graphics.off()
pdf("graphics/undercoverage_lrv_15.pdf",pointsize=16,width=8.0,height=6.0)
plot(NULL, xlab="AR(1) coefficient", ylab="Coverage", main=expression(paste(T[0],"=15")), col="blue", pch=".", xlim = c(0,1), ylim=c(0.4,1))
lines(rho.vec,results_overall_t_15,col="black",lwd=5, lty=1,type="b",pch=1)
lines(rho.vec,results_overall_nw_15,col="gray",lwd=5, lty=1,type="b",pch=2)
abline(h=0.9,lty=1,lwd=1,col="gray")
legend("bottomleft",legend=c("t-Test (K=3)","Conventional Newey-West standard errors"), seg.len=2, col=c("black","gray"),fill=NA,border=NA, lty=c(1,1),lwd=c(5,5), pch=c(1,2), merge=T,bty="n")
dev.off()

graphics.off()
pdf("graphics/undercoverage_lrv_30.pdf",pointsize=16,width=8.0,height=6.0)
plot(NULL, xlab="AR(1) coefficient", ylab="Coverage", main=expression(paste(T[0],"=30")), col="blue", pch=".", xlim = c(0,1), ylim=c(0.4,1))
lines(rho.vec,results_overall_t_30,col="black",lwd=5, lty=1,type="b",pch=1)
lines(rho.vec,results_overall_nw_30,col="gray",lwd=5, lty=1,type="b",pch=2)
abline(h=0.9,lty=1,lwd=1,col="gray")
legend("bottomleft",legend=c("t-Test (K=3)","Conventional Newey-West standard errors"), seg.len=2, col=c("black","gray"),fill=NA,border=NA, lty=c(1,1),lwd=c(5,5), pch=c(1,2), merge=T,bty="n")
dev.off()

